/*
 * mmd.c
 *
 * **************************************************************
 * The following C function was developed from a FORTRAN subroutine
 * in SPARSPAK written by Eleanor Chu, Alan George, Joseph Liu
 * and Esmond Ng.
 * 
 * The FORTRAN-to-C transformation and modifications such as dynamic
 * memory allocation and deallocation were performed by Chunguang
 * Sun.
 * ************************************************************** 
 *
 * Taken from SMMS, George 12/13/94
 *
 * The meaning of invperm, and perm vectors is different from that
 * in genqmd_ of SparsPak
 *
 * $Id: mmd.c,v 1.1 1998/11/27 17:59:25 karypis Exp $
 */

namespace cusp
{
namespace graph
{

/*************************************************************************
*  genmmd  -- multiple minimum external degree
*  purpose -- this routine implements the minimum degree
*     algorithm. it makes use of the implicit representation
*     of elimination graphs by quotient graphs, and the notion
*     of indistinguishable nodes. It also implements the modifications
*     by multiple elimination and minimum external degree.
*     Caution -- the adjacency vector adjncy will be destroyed.
*  Input parameters --
*     neqns -- number of equations.
*     (xadj, adjncy) -- the adjacency structure.
*     delta  -- tolerance value for multiple elimination.
*     maxint -- maximum machine representable (short) integer
*               (any smaller estimate will do) for marking nodes.
*  Output parameters --
*     perm -- the minimum degree ordering.
*     invp -- the inverse of perm.
*     *ncsub -- an upper bound on the number of nonzero subscripts
*               for the compressed storage scheme.
*  Working parameters --
*     head -- vector for head of degree lists.
*     invp  -- used temporarily for degree forward link.
*     perm  -- used temporarily for degree backward link.
*     qsize -- vector for size of supernodes.
*     list -- vector for temporary linked lists.
*     marker -- a temporary marker vector.
*  Subroutines used -- mmdelm, mmdint, mmdnum, mmdupd.
**************************************************************************/
void genmmd(int neqns, int *xadj, int *adjncy, int *invp, int *perm,
	    int delta, int *head, int *qsize, int *list, int *marker,
	    int maxint, int *ncsub)
{
  int  ehead, i, mdeg, mdlmt, mdeg_node, nextmd, num, tag;
/* ------------------------------prototypes---------------------------- */ 
  int  mmdint(int, int*, int*, int*,int*, int*, int*, int*, int*);
  void mmdelm(int, int*, int*, int*, int*, int*, int*, int*, int*,int,int);
  void mmdupd(int, int, int*, int*, int, int*, int*, int*, int*,
	      int*, int*, int*, int,int*); 
  void mmdnum(int, int*, int*, int*); 
/*-----------------------------------------------------------------------*/

  if (neqns <= 0)  
    return;
  
  /* Adjust from C to Fortran */
  xadj--; adjncy--; invp--; perm--; head--; qsize--; list--; marker--;
  
  /* initialization for the minimum degree algorithm. */
  *ncsub = 0;
  mmdint(neqns, xadj, adjncy, head, invp, perm, qsize, list, marker);
  
  /*  'num' counts the number of ordered nodes plus 1. */
  num = 1;
  
  /* eliminate all isolated nodes. */
  nextmd = head[1];
  while (nextmd > 0) {
    mdeg_node = nextmd;
    nextmd = invp[mdeg_node];
    marker[mdeg_node] = maxint;
    invp[mdeg_node] = -num;
    num = num + 1;
  }
  
  /* search for node of the minimum degree. 'mdeg' is the current */
  /* minimum degree; 'tag' is used to facilitate marking nodes.   */
  if (num > neqns) 
    goto n1000;
  tag = 1;
  head[1] = 0;
  mdeg = 2;
  
  /* infinite loop here ! */
  while (1) {
    while (head[mdeg] <= 0) 
      mdeg++;
    
    /* use value of 'delta' to set up 'mdlmt', which governs */
    /* when a degree update is to be performed.              */
    mdlmt = mdeg + delta;
    ehead = 0;
    
  n500:
    mdeg_node = head[mdeg];
    while (mdeg_node <= 0) {
      mdeg++;
      
      if (mdeg > mdlmt) 
	goto n900;
      mdeg_node = head[mdeg];
    };
    
    /*  remove 'mdeg_node' from the degree structure. */
    nextmd = invp[mdeg_node];
    head[mdeg] = nextmd;
    if (nextmd > 0)  
      perm[nextmd] = -mdeg;
    invp[mdeg_node] = -num;
    *ncsub += mdeg + qsize[mdeg_node] - 2;
    if ((num+qsize[mdeg_node]) > neqns)  
      goto n1000;
    
    /*  eliminate 'mdeg_node' and perform quotient graph */
    /*  transformation. reset 'tag' value if necessary.    */
    tag++;
    if (tag >= maxint) {
      tag = 1;
      for (i = 1; i <= neqns; i++)
	if (marker[i] < maxint)  
	  marker[i] = 0;
    };
    
    mmdelm(mdeg_node, xadj, adjncy, head, invp, perm,
	   qsize, list, marker, maxint, tag);
    
    num += qsize[mdeg_node];
    list[mdeg_node] = ehead;
    ehead = mdeg_node;
    if (delta >= 0) 
      goto n500;
    
  n900:
    /* update degrees of the nodes involved in the  */
    /* minimum degree nodes elimination.            */
    if (num > neqns)  
      goto n1000;
    mmdupd( ehead, neqns, xadj, adjncy, delta, &mdeg, head, invp, perm, qsize, list, marker, maxint, &tag);
  }; /* end of -- while ( 1 ) -- */
  
 n1000:
  mmdnum( neqns, perm, invp, qsize );
  
  /* Adjust from Fortran back to C*/
  xadj++; adjncy++; invp++; perm++; head++; qsize++; list++; marker++;
}

/**************************************************************************
*           mmdelm ...... multiple minimum degree elimination
* Purpose -- This routine eliminates the node mdeg_node of minimum degree
*     from the adjacency structure, which is stored in the quotient
*     graph format. It also transforms the quotient graph representation
*     of the elimination graph.
* Input parameters --
*     mdeg_node -- node of minimum degree.
*     maxint -- estimate of maximum representable (short) integer.
*     tag    -- tag value.
* Updated parameters --
*     (xadj, adjncy) -- updated adjacency structure.
*     (head, forward, backward) -- degree doubly linked structure.
*     qsize -- size of supernode.
*     marker -- marker vector.
*     list -- temporary linked list of eliminated nabors.
***************************************************************************/
void mmdelm(int mdeg_node, int *xadj, int *adjncy, int *head, int *forward,
	    int *backward, int *qsize, int *list, int *marker, int maxint,
	    int tag)
{
  int   element, i,   istop, istart, j,
    jstop, jstart, link,
    nabor, node, npv, nqnbrs, nxnode,
    pvnode, rlmt, rloc, rnode, xqnbr;
  
  /* find the reachable set of 'mdeg_node' and */
  /* place it in the data structure.           */
  marker[mdeg_node] = tag;
  istart = xadj[mdeg_node];
  istop = xadj[mdeg_node+1] - 1;
  
  /* 'element' points to the beginning of the list of  */
  /* eliminated nabors of 'mdeg_node', and 'rloc' gives the */
  /* storage location for the next reachable node.   */
  element = 0;
  rloc = istart;
  rlmt = istop;
  for ( i = istart; i <= istop; i++ ) {
    nabor = adjncy[i];
    if ( nabor == 0 ) break;
    if ( marker[nabor] < tag ) {
      marker[nabor] = tag;
      if ( forward[nabor] < 0 )  {
	list[nabor] = element;
	element = nabor;
      } else {
	adjncy[rloc] = nabor;
	rloc++;
      };
    }; /* end of -- if -- */
  }; /* end of -- for -- */
  
  /* merge with reachable nodes from generalized elements. */
  while ( element > 0 ) {
    adjncy[rlmt] = -element;
    link = element;
    
  n400:
    jstart = xadj[link];
    jstop = xadj[link+1] - 1;
    for ( j = jstart; j <= jstop; j++ ) {
      node = adjncy[j];
      link = -node;
      if ( node < 0 )  goto n400;
      if ( node == 0 ) break;
      if ((marker[node]<tag)&&(forward[node]>=0)) {
	marker[node] = tag;
	/*use storage from eliminated nodes if necessary.*/
	while ( rloc >= rlmt ) {
	  link = -adjncy[rlmt];
	  rloc = xadj[link];
	  rlmt = xadj[link+1] - 1;
	};
	adjncy[rloc] = node;
	rloc++;
      };
    }; /* end of -- for ( j = jstart; -- */
    element = list[element];
  };  /* end of -- while ( element > 0 ) -- */
  if ( rloc <= rlmt ) adjncy[rloc] = 0;
  /* for each node in the reachable set, do the following. */
  link = mdeg_node;
  
 n1100:
  istart = xadj[link];
  istop = xadj[link+1] - 1;
  for ( i = istart; i <= istop; i++ ) {
    rnode = adjncy[i];
    link = -rnode;
    if ( rnode < 0 ) goto n1100;
    if ( rnode == 0 ) return;
    
    /* 'rnode' is in the degree list structure. */
    pvnode = backward[rnode];
    if (( pvnode != 0 ) && ( pvnode != (-maxint) )) {
      /* then remove 'rnode' from the structure. */
      nxnode = forward[rnode];
      if ( nxnode > 0 ) backward[nxnode] = pvnode;
      if ( pvnode > 0 ) forward[pvnode] = nxnode;
      npv = -pvnode;
      if ( pvnode < 0 ) head[npv] = nxnode;
    };
    
    /* purge inactive quotient nabors of 'rnode'. */
    jstart = xadj[rnode];
    jstop = xadj[rnode+1] - 1;
    xqnbr = jstart;
    for ( j = jstart; j <= jstop; j++ ) {
      nabor = adjncy[j];
      if ( nabor == 0 ) break;
      if ( marker[nabor] < tag ) {
	adjncy[xqnbr] = nabor;
	xqnbr++;
      };
    };
    
    /* no active nabor after the purging. */
    nqnbrs = xqnbr - jstart;
    if ( nqnbrs <= 0 ) {
      /* merge 'rnode' with 'mdeg_node'. */
      qsize[mdeg_node] += qsize[rnode];
      qsize[rnode] = 0;
      marker[rnode] = maxint;
      forward[rnode] = -mdeg_node;
      backward[rnode] = -maxint;
    } else {
      /* flag 'rnode' for degree update, and  */
      /* add 'mdeg_node' as a nabor of 'rnode'.      */
      forward[rnode] = nqnbrs + 1;
      backward[rnode] = 0;
      adjncy[xqnbr] = mdeg_node;
      xqnbr++;
      if ( xqnbr <= jstop )  adjncy[xqnbr] = 0;
    };
  }; /* end of -- for ( i = istart; -- */
  return;
}

/***************************************************************************
*    mmdint ---- mult minimum degree initialization
*    purpose -- this routine performs initialization for the
*       multiple elimination version of the minimum degree algorithm.
*    input parameters --
*       neqns  -- number of equations.
*       (xadj, adjncy) -- adjacency structure.
*    output parameters --
*       (head, dfrow, backward) -- degree doubly linked structure.
*       qsize -- size of supernode ( initialized to one).
*       list -- linked list.
*       marker -- marker vector.
****************************************************************************/
int  mmdint(int neqns, int *xadj, int *adjncy, int *head, int *forward,
	    int *backward, int *qsize, int *list, int *marker)
{
  int  fnode, ndeg, node;
  
  for ( node = 1; node <= neqns; node++ ) {
    head[node] = 0;
    qsize[node] = 1;
    marker[node] = 0;
    list[node] = 0;
  };
  
  /* initialize the degree doubly linked lists. */
  for ( node = 1; node <= neqns; node++ ) {
    ndeg = xadj[node+1] - xadj[node]/* + 1*/;   /* george */
    if (ndeg == 0)
      ndeg = 1;
    fnode = head[ndeg];
    forward[node] = fnode;
    head[ndeg] = node;
    if ( fnode > 0 ) backward[fnode] = node;
    backward[node] = -ndeg;
  };
  return 0;
}

/****************************************************************************
 * mmdnum --- multi minimum degree numbering
 * purpose -- this routine performs the final step in producing
 *    the permutation and inverse permutation vectors in the
 *    multiple elimination version of the minimum degree
 *    ordering algorithm.
 * input parameters --
 *     neqns -- number of equations.
 *     qsize -- size of supernodes at elimination.
 * updated parameters --
 *     invp -- inverse permutation vector. on input,
 *             if qsize[node] = 0, then node has been merged
 *             into the node -invp[node]; otherwise,
 *            -invp[node] is its inverse labelling.
 * output parameters --
 *     perm -- the permutation vector.
 ****************************************************************************/
void mmdnum(int neqns, int *perm, int *invp, int *qsize)
{
  int father, nextf, node, nqsize, num, root;
  
  for ( node = 1; node <= neqns; node++ ) {
    nqsize = qsize[node];
    if ( nqsize <= 0 ) perm[node] = invp[node];
    if ( nqsize > 0 )  perm[node] = -invp[node];
  };
  
  /* for each node which has been merged, do the following. */
  for ( node = 1; node <= neqns; node++ ) {
    if ( perm[node] <= 0 )  {
      
      /* trace the merged tree until one which has not */
      /* been merged, call it root.                    */
      father = node;
      while ( perm[father] <= 0 )
	father = - perm[father];
      
      /* number node after root. */
      root = father;
      num = perm[root] + 1;
      invp[node] = -num;
      perm[root] = num;

      /* shorten the merged tree. */
      father = node;
      nextf = - perm[father];
      while ( nextf > 0 ) {
	perm[father] = -root;
	father = nextf;
	nextf = -perm[father];
      };
    };  /* end of -- if ( perm[node] <= 0 ) -- */
  }; /* end of -- for ( node = 1; -- */
  
  /* ready to compute perm. */
  for ( node = 1; node <= neqns; node++ ) {
    num = -invp[node];
    invp[node] = num;
    perm[num] = node;
  };
  return;
}

/****************************************************************************
* mmdupd ---- multiple minimum degree update
* purpose -- this routine updates the degrees of nodes after a
*            multiple elimination step.
* input parameters --
*    ehead -- the beginning of the list of eliminated nodes
*             (i.e., newly formed elements).
*    neqns -- number of equations.
*    (xadj, adjncy) -- adjacency structure.
*    delta -- tolerance value for multiple elimination.
*    maxint -- maximum machine representable (short) integer.
* updated parameters --
*    mdeg -- new minimum degree after degree update.
*    (head, forward, backward) -- degree doubly linked structure.
*    qsize -- size of supernode.
*    list -- marker vector for degree update.
*    *tag   -- tag value.
****************************************************************************/
void mmdupd(int ehead, int neqns, int *xadj, int *adjncy, int delta, int *mdeg,
	    int *head, int *forward, int *backward, int *qsize, int *list,
	    int *marker, int maxint,int *tag)
{
  int  deg, deg0, element, enode, fnode, i, iq2, istop,
    istart, j, jstop, jstart, link, mdeg0, mtag, nabor,
    node, q2head, qxhead;
  
  mdeg0 = *mdeg + delta;
  element = ehead;
  
 n100:
  if ( element <= 0 ) return;
  
  /* for each of the newly formed element, do the following. */
  /* reset tag value if necessary.                           */
  mtag = *tag + mdeg0;
  if ( mtag >= maxint ) {
    *tag = 1;
    for ( i = 1; i <= neqns; i++ )
      if ( marker[i] < maxint ) marker[i] = 0;
    mtag = *tag + mdeg0;
  };
  
  /* create two linked lists from nodes associated with 'element': */
  /* one with two nabors (q2head) in the adjacency structure, and the*/
  /* other with more than two nabors (qxhead). also compute 'deg0',*/
  /* number of nodes in this element.                              */
  q2head = 0;
  qxhead = 0;
  deg0 = 0;
  link =element;
  
 n400:
  istart = xadj[link];
  istop = xadj[link+1] - 1;
  for ( i = istart; i <= istop; i++ ) {
    enode = adjncy[i];
    link = -enode;
    if ( enode < 0 )  goto n400;
    if ( enode == 0 ) break;
    if ( qsize[enode] != 0 ) {
      deg0 += qsize[enode];
      marker[enode] = mtag;
      
      /*'enode' requires a degree update*/
      if ( backward[enode] == 0 ) {
	/* place either in qxhead or q2head list. */
	if ( forward[enode] != 2 ) {
	  list[enode] = qxhead;
	  qxhead = enode;
	} else {
	  list[enode] = q2head;
	  q2head = enode;
	};
      };
    }; /* enf of -- if ( qsize[enode] != 0 ) -- */
  }; /* end of -- for ( i = istart; -- */
  
  /* for each node in q2 list, do the following. */
  enode = q2head;
  iq2 = 1;
  
 n900:
  if ( enode <= 0 ) goto n1500;
  if ( backward[enode] != 0 ) goto n2200;
  (*tag)++;
  deg = deg0;
  
  /* identify the other adjacent element nabor. */
  istart = xadj[enode];
  nabor = adjncy[istart];
  if ( nabor == element ) nabor = adjncy[istart+1];
  link = nabor;
  if ( forward[nabor] >= 0 ) {
    /* nabor is uneliminated, increase degree count. */
    deg += qsize[nabor];
    goto n2100;
  };
  
  /* the nabor is eliminated. for each node in the 2nd element */
  /* do the following.                                         */
 n1000:
  istart = xadj[link];
  istop = xadj[link+1] - 1;
       for ( i = istart; i <= istop; i++ ) {
	 node = adjncy[i];
	 link = -node;
	 if ( node != enode ) {
	   if ( node < 0 ) goto n1000;
	   if ( node == 0 )  goto n2100;
	   if ( qsize[node] != 0 ) {
	     if ( marker[node] < *tag ) {
	       /* 'node' is not yet considered. */
	       marker[node] = *tag;
	       deg += qsize[node];
	     } else {
	       if ( backward[node] == 0 ) {
		 if ( forward[node] == 2 ) {
                                /* 'node' is indistinguishable from 'enode'.*/
                                /* merge them into a new supernode.         */
		   qsize[enode] += qsize[node];
		   qsize[node] = 0;
		   marker[node] = maxint;
		   forward[node] = -enode;
		   backward[node] = -maxint;
		 } else {
                                /* 'node' is outmacthed by 'enode' */
		   if (backward[node]==0) backward[node] = -maxint;
		 };
	       }; /* end of -- if ( backward[node] == 0 ) -- */
	     }; /* end of -- if ( marker[node] < *tag ) -- */
	   }; /* end of -- if ( qsize[node] != 0 ) -- */
	 }; /* end of -- if ( node != enode ) -- */
       }; /* end of -- for ( i = istart; -- */
       goto n2100;
       
 n1500:
       /* for each 'enode' in the 'qx' list, do the following. */
       enode = qxhead;
       iq2 = 0;
       
 n1600:    if ( enode <= 0 )  goto n2300;
       if ( backward[enode] != 0 )  goto n2200;
       (*tag)++;
       deg = deg0;

       /*for each unmarked nabor of 'enode', do the following.*/
       istart = xadj[enode];
       istop = xadj[enode+1] - 1;
       for ( i = istart; i <= istop; i++ ) {
	 nabor = adjncy[i];
	 if ( nabor == 0 ) break;
	 if ( marker[nabor] < *tag ) {
	   marker[nabor] = *tag;
	   link = nabor;
	   if ( forward[nabor] >= 0 ) 
	     /*if uneliminated, include it in deg count.*/
	     deg += qsize[nabor];
	   else {
n1700:
	     /* if eliminated, include unmarked nodes in this*/
	     /* element into the degree count.             */
	     jstart = xadj[link];
	     jstop = xadj[link+1] - 1;
	     for ( j = jstart; j <= jstop; j++ ) {
	       node = adjncy[j];
	       link = -node;
	       if ( node < 0 ) goto n1700;
	       if ( node == 0 ) break;
	       if ( marker[node] < *tag ) {
		 marker[node] = *tag;
		 deg += qsize[node];
	       };
	     }; /* end of -- for ( j = jstart; -- */
	   }; /* end of -- if ( forward[nabor] >= 0 ) -- */
	 }; /* end of -- if ( marker[nabor] < *tag ) -- */
       }; /* end of -- for ( i = istart; -- */
       
 n2100:
       /* update external degree of 'enode' in degree structure, */
       /* and '*mdeg' if necessary.                     */
       deg = deg - qsize[enode] + 1;
       fnode = head[deg];
       forward[enode] = fnode;
       backward[enode] = -deg;
       if ( fnode > 0 ) backward[fnode] = enode;
       head[deg] = enode;
       if ( deg < *mdeg ) *mdeg = deg;
       
 n2200:
       /* get next enode in current element. */
       enode = list[enode];
       if ( iq2 == 1 ) goto n900;
       goto n1600;
       
 n2300:
       /* get next element in the list. */
       *tag = mtag;
       element = list[element];
       goto n100;
}

} // end namespace graph
} // end namespace cusp
